Introduction

Terrestrial Ecology Department, Technical University of Munich

In a greenhouse full factorial experiment, we investigated the effects of prothioconazole and fluoxastrobin fungicide irrigation, of inoculation of the pathogenic fungus Drechslera teres and of barley cultivars on plant growth and grain yield variables, and particularly on barley root endophytic communities. Our experimental design was an approach to real agriculture practice, where generally fungicides are irrigated prophylactic before occurrence of pathogenic fungi. Additionally, aphid abundance has been measured. 60 root samples with all treatment combinations were chosen for DNA-extraction, PCR amplification and metabarcoding. For beta diversity analysis, multivariate statistics had been performed to reduce the large data sets into few principal coordinates in a MDS-CAP ordination. We hypothesized that the fungicide treatment will affect alpha and beta diversity in bacterial and fungal endophyte communities. It was found, that the fungicide and Drechslera treatments had no significant effect on alpha and beta diversity of the endophytes, but barley cultivar showed to be a significant driver of beta diversity of bacterial and fungal barley root endophyte communities.

Bacteria

Setting envriroment

loading packages

##Loading packages

library(phyloseq)
library(ggplot2)    
library(readxl)
library(reshape2)
library(dplyr)        
library(microbiome)
library(DESeq2)
library(writexl)
library(microeco)
library(vegan)
library(microbiomeutilities)
library(lme4)
library(lmerTest)
library(nlme)
library(metagMisc)
library(gridExtra)
library(microeco)
library(file2meco)
library(ggalluvial)
library(ggcorrplot)
library(ggpubr)
library(car)
library(DT)


## Setting theme 

theme_set(theme_classic())

Loading data

setwd("C:/Users/yurip/Desktop/COMI/Colaborations/Stephan/second_round")

otu_mat<- as.data.frame(read_excel("asvs.xlsx"))
tax_mat<-  
  as.data.frame(read_excel("taxonomy.xlsx",
                           na = "NA",
                           col_names =
                             c("ASV", "Kingdom", "Phylum", "Class", "Order", "Family", "Genus", 
                                         "Species"))) %>%
  purrr::map_df(~gsub(
    pattern = "metagenome|uncultured|unidentified|Unknown",
    replacement = NA,
    .x)) %>%
  mutate_if(is.character, stringr::str_trim) %>%
  mutate(Kingdom = ifelse(is.na(Kingdom),
                          "U. Kingdom",
                          Kingdom),
         Phylum = coalesce(Phylum,
                           ifelse(grepl("^U.", Kingdom),
                                  Kingdom,
                                  paste("K_", Kingdom))),
         Class = coalesce(Class,
                          ifelse(grepl("^U.", Phylum),
                                 Phylum,
                                 paste("P_", Phylum))),
         Order = coalesce(Order,
                          ifelse(grepl("^U.", Class),
                                 Class,
                                 paste("C_", Class))),
         Family = coalesce(Family,
                           ifelse(grepl("^U.", Order),
                                  Order,
                                  paste("O_", Order))),
         Genus = coalesce(Genus,
                          ifelse(grepl("^U.", Family),
                                 Family,
                                 paste("F_", Family))),
         Species = coalesce(Species,
                            ifelse(grepl("^U.", Genus),
                                   Genus,
                                   paste("G_", Genus)))) %>%
  tibble::column_to_rownames("ASV") %>%
  filter(Kingdom %in% "Bacteria" &
           !grepl("Chloroplast", Genus) &
           !grepl("Mitochondria", Genus) &
           !Order %in% "Chloroplast")

samples_df <- as.data.frame(read_excel("metadata_stephan2.xlsx"))

otu_names=otu_mat$OTUS
row.names(otu_mat) <- otu_names
otu_mat <- otu_mat[,2:61]


row.names(samples_df) <- samples_df$sample
samples_df <- samples_df %>% select (-sample) 
samples_df$Block <- factor(samples_df$Block, levels = c("Block1", "Block2", "Block3", "Block4","Block5" ))


otu_mat <- as.matrix(otu_mat)
tax_mat <- as.matrix(tax_mat)

class(otu_mat) <- "numeric"

OTU = otu_table(otu_mat, taxa_are_rows = TRUE)
TAX = tax_table(tax_mat)
samples = sample_data(samples_df)

stephan <- phyloseq(OTU, TAX, samples)
stephan
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 6275 taxa and 60 samples ]
## sample_data() Sample Data:       [ 60 samples by 11 sample variables ]
## tax_table()   Taxonomy Table:    [ 6275 taxa by 7 taxonomic ranks ]
sample_names(stephan)
##  [1] "Weisser.1_S34"   "Weisser.3_S35"   "Weisser.5_S36"   "Weisser.7_S37"  
##  [5] "Weisser.9_S38"   "Weisser.11_S39"  "Weisser.13_S40"  "Weisser.15_S41" 
##  [9] "Weisser.17_S42"  "Weisser.19_S43"  "Weisser.21_S44"  "Weisser.23_S45" 
## [13] "Weisser.25_S46"  "Weisser.27_S47"  "Weisser.29_S48"  "Weisser.31_S49" 
## [17] "Weisser.33_S50"  "Weisser.35_S51"  "Weisser.37_S52"  "Weisser.39_S53" 
## [21] "Weisser.41_S54"  "Weisser.43_S55"  "Weisser.45_S56"  "Weisser.47_S57" 
## [25] "Weisser.49_S58"  "Weisser.51_S59"  "Weisser.53_S60"  "Weisser.55_S61" 
## [29] "Weisser.57_S62"  "Weisser.59_S63"  "Weisser.61_S64"  "Weisser.63_S65" 
## [33] "Weisser.65_S66"  "Weisser.67_S67"  "Weisser.69_S68"  "Weisser.71_S69" 
## [37] "Weisser.73_S70"  "Weisser.75_S71"  "Weisser.77_S72"  "Weisser.79_S73" 
## [41] "Weisser.81_S74"  "Weisser.83_S75"  "Weisser.85_S76"  "Weisser.87_S77" 
## [45] "Weisser.89_S78"  "Weisser.91_S79"  "Weisser.93_S80"  "Weisser.95_S81" 
## [49] "Weisser.97_S82"  "Weisser.99_S83"  "Weisser.101_S84" "Weisser.103_S85"
## [53] "Weisser.105_S86" "Weisser.107_S87" "Weisser.109_S88" "Weisser.111_S89"
## [57] "Weisser.113_S90" "Weisser.115_S91" "Weisser.117_S92" "Weisser.119_S93"
rank_names(stephan)
## [1] "Kingdom" "Phylum"  "Class"   "Order"   "Family"  "Genus"   "Species"
sample_variables(stephan)
##  [1] "Plant_number"             "Cultivar"                
##  [3] "Fungicide_treatment"      "Drechslera_treatment"    
##  [5] "Aphid_treatment"          "Block"                   
##  [7] "Cultivar_and_Drechslera"  "Fungicide_and_Drechslera"
##  [9] "Aphid_count"              "Total_Seed_Biomass"      
## [11] "Dry_Weight"
rank_names(stephan)
## [1] "Kingdom" "Phylum"  "Class"   "Order"   "Family"  "Genus"   "Species"
sample_variables(stephan)
##  [1] "Plant_number"             "Cultivar"                
##  [3] "Fungicide_treatment"      "Drechslera_treatment"    
##  [5] "Aphid_treatment"          "Block"                   
##  [7] "Cultivar_and_Drechslera"  "Fungicide_and_Drechslera"
##  [9] "Aphid_count"              "Total_Seed_Biomass"      
## [11] "Dry_Weight"

Rarefaction curve and normalization

mat <- t(otu_table(stephan))
class(mat) <- "matrix"
rarecurve(mat, step=50, cex=0.5, label = TRUE, 
          xlab = "Reads", ylab = "ASVs")

#Normalize number of reads in each sample using median sequencing depth.
total = median(sample_sums(stephan))
standf = function(x, t=total) round(t * (x / sum(x)))
stephan.rarefied = transform_sample_counts(stephan, standf)

Diversity analyzes

Alpha diversity

alpha_bac=plot_richness(stephan.rarefied, x="Cultivar", color="Fungicide_treatment",
                        measures=c("Shannon")) + geom_boxplot() +
  facet_wrap(~Drechslera_treatment, scales = "free_x", nrow = 1) +
  scale_color_manual(values = c("#90EE90", "#228822")) +
  labs(title = "Alpha diversity (Shannon) - Bacteria")

alpha_bac

alpha_bac2=plot_richness(stephan.rarefied, x="Cultivar", color="Fungicide_treatment",
                        measures=c("Simpson")) + geom_boxplot() +
  facet_wrap(~Drechslera_treatment, scales = "free_x", nrow = 1) +
  scale_color_manual(values = c("#90EE90", "#228822")) +
  labs(title = "Alpha diversity (Simpson) - Bacteria")

alpha_bac2

## Calculation diversity metrics 

diversity=estimate_richness(stephan.rarefied)
shannon=diversity$Shannon
observed=diversity$Observed
observed.log <- log(observed)
evenness=shannon/observed.log
diversity$evennes=evenness
richness=as.data.frame(diversity)


DT::datatable(richness, class = "cell-border stripe", rownames = T,
              filter = "top", editable = TRUE, extensions = 'Buttons', options = list(
  dom = 'Bfrtip',
  buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
  scrollY = '300px',
  scrollCollapse = TRUE
))

Beta diversity

### PCoA and Constrained Analysis of Principal Coordinates (CAP)

ps1V4.rel_bray <- phyloseq::distance(stephan.rarefied, method = "bray") # CAP ordinate 
cap_ord <- ordinate(physeq = stephan.rarefied,  
                    method = "CAP", 
                    distance = ps1V4.rel_bray, 
                    formula = ~ Aphid_count + Total_Seed_Biomass +
                      Dry_Weight)

scree.cap <- plot_scree(cap_ord, "Scree Plot for MCs in Constrained Analysis of Principal Coordinates - CAPSCALE")
print(scree.cap)

cap_plot <- plot_ordination(physeq = stephan.rarefied, 
                            ordination = cap_ord, 
                            color = "Cultivar", 
                            shape = "Fungicide_treatment",
                            axes = c(1,2)) + 
  geom_point(aes(colour = Cultivar), size = 3) + 
  geom_point(size = 3) +
  scale_color_manual(
    values = c( "#90EE90", "#228822", "#5F7FC7", "orange","#DA5724")) 



#Now add the environmental variables as arrows
cap_ord
## Call: capscale(formula = distance ~ Aphid_count + Total_Seed_Biomass +
## Dry_Weight, data = data)
## 
##                  Inertia Proportion Rank
## Total         18.9706062  1.0000000     
## Constrained    1.3895020  0.0732450    3
## Unconstrained 17.5853700  0.9269799   56
## Imaginary     -0.0042658 -0.0002249    1
## Inertia is squared Bray distance 
## 
## Eigenvalues for constrained axes:
##   CAP1   CAP2   CAP3 
## 0.6710 0.4388 0.2797 
## 
## Eigenvalues for unconstrained axes:
##  MDS1  MDS2  MDS3  MDS4  MDS5  MDS6  MDS7  MDS8 
## 3.194 1.321 1.047 0.874 0.732 0.663 0.591 0.536 
## (Showing 8 of 56 unconstrained eigenvalues)
arrowmat <- vegan::scores(cap_ord, display = "bp")

# Add labels, make a data.frame
arrowdf <- data.frame(labels = rownames(arrowmat), arrowmat)

# Define the arrow aesthetic mapping
arrow_map <- aes(xend = CAP1, 
                 yend = CAP1, 
                 x = 0, 
                 y = 0, 
                 shape = NULL, 
                 color = NULL, 
                 label = labels)

label_map <- aes(x = 1.3 * CAP1, 
                 y = 1.3 * CAP1, 
                 shape = NULL, 
                 color = NULL, 
                 label = labels)

arrowhead = arrow(length = unit(0.02, "npc"))

# Make a new graphic
cap_plot + 
  geom_segment(
    mapping = arrow_map, 
    size = .5, 
    data = arrowdf, 
    color = "black", 
    arrow = arrowhead
  ) + 
  geom_text(
    mapping = label_map, 
    size = 4,  
    data = arrowdf, 
    show.legend = FALSE
  ) +
  stat_ellipse(type = "t") + ggtitle("CAP_Plot - Bacteria")

Taxonomic composition

Box plot

mycols <- c("yellow", "#90EE90", "#228822", "#5F7FC7", "orange","#DA5724")

pn <- plot_taxa_boxplot(stephan.rarefied,
                        taxonomic.level = "Family",
                        top.otu = 9,
                        group = "Cultivar_and_Drechslera",
                        title = "top 9 taxa at family level - Bacteria",
                        keep.other = FALSE,
                        add.violin= FALSE,
                        group.order = c("ACCO Drechs_NO", "ACCO Drechs_YES",
                                        "KLAR Drechs_NO", "KLAR Drechs_YES" ,
                                        "PLAN Drechs_YES", "PLAN Drechs_NO"),
                        group.colors = mycols
)
   
box_plot_bacteria = pn + theme_biome_utils() +  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  theme(legend.position = "bottom")

box_plot_bacteria

### Exporting data most abundant taxa abundances for linear models ###

top_nine_bacteria <- subset_taxa(stephan.rarefied, Family %in% c("Cellvibrionaceae" ,  "Comamonadaceae",   
                                                                 "Enterobacteriaceae", "Oxalobacteraceae" , 
                                                                 "Pseudomonadaceae", "Rhizobiaceae",       
                                                                 "Rhodanobacteraceae", "Xanthobacteraceae",  
                                                                 "Xanthomonadaceae"))

top_nine_bacteria
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1298 taxa and 60 samples ]
## sample_data() Sample Data:       [ 60 samples by 11 sample variables ]
## tax_table()   Taxonomy Table:    [ 1298 taxa by 7 taxonomic ranks ]
top_nine_melted=psmelt(top_nine_bacteria)


write_xlsx(top_nine_melted,"top_nine_bacteria.xlsx")

Heatmap

meco_dataset <- phyloseq2meco(stephan.rarefied)

t1 <- trans_abund$new(dataset = meco_dataset, taxrank = "Family", ntaxa = 10, groupmean = "Fungicide_and_Drechslera")

t1 <- trans_abund$new(dataset = meco_dataset, taxrank = "Family", ntaxa = 20)
t1$plot_heatmap(facet = "Cultivar", xtext_keep = FALSE, withmargin = FALSE)

Statistical analyzes

Permanova

stephan_bray <- phyloseq::distance(stephan, method = "bray")
stephandf <- data.frame(sample_data(stephan))

adonis2(stephan_bray ~ Fungicide_treatment*Drechslera_treatment*Cultivar, data = stephandf)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = stephan_bray ~ Fungicide_treatment * Drechslera_treatment * Cultivar, data = stephandf)
##                                                   Df SumOfSqs      R2      F
## Fungicide_treatment                                1   0.3135 0.01637 0.9617
## Drechslera_treatment                               1   0.2633 0.01375 0.8077
## Cultivar                                           2   1.0527 0.05498 1.6146
## Fungicide_treatment:Drechslera_treatment           1   0.2451 0.01280 0.7519
## Fungicide_treatment:Cultivar                       2   0.4920 0.02570 0.7546
## Drechslera_treatment:Cultivar                      2   0.5335 0.02786 0.8183
## Fungicide_treatment:Drechslera_treatment:Cultivar  2   0.5986 0.03126 0.9181
## Residual                                          48  15.6482 0.81727       
## Total                                             59  19.1470 1.00000       
##                                                   Pr(>F)   
## Fungicide_treatment                                0.464   
## Drechslera_treatment                               0.754   
## Cultivar                                           0.009 **
## Fungicide_treatment:Drechslera_treatment           0.858   
## Fungicide_treatment:Cultivar                       0.929   
## Drechslera_treatment:Cultivar                      0.830   
## Fungicide_treatment:Drechslera_treatment:Cultivar  0.624   
## Residual                                                   
## Total                                                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Linear mixed effect models

richness$Block = samples_df$Block
richness$Drechslera_treatment = samples_df$Drechslera_treatment
richness$Fungicide_treatment = samples_df$Fungicide_treatment
richness$Sample = samples_df$Plant_number
richness$Cultivar = samples_df$Cultivar


table(richness$Fungicide_treatment,richness$Drechslera_treatment)
##                
##                 Drechs_NO Drechs_YES
##   Fungicide_NO         15         15
##   Fungicide_YES        15         15
m1<-lme(Shannon~Drechslera_treatment*Fungicide_treatment*Cultivar,                            ,
        random=~1|Block,
        na.action=na.exclude,
        data=richness)

anova(m1)
##                                                   numDF denDF   F-value p-value
## (Intercept)                                           1    44 1158.7360  <.0001
## Drechslera_treatment                                  1    44    1.6947  0.1998
## Fungicide_treatment                                   1    44    0.0191  0.8908
## Cultivar                                              2    44    0.8794  0.4222
## Drechslera_treatment:Fungicide_treatment              1    44    0.5659  0.4559
## Drechslera_treatment:Cultivar                         2    44    0.9959  0.3776
## Fungicide_treatment:Cultivar                          2    44    0.8806  0.4217
## Drechslera_treatment:Fungicide_treatment:Cultivar     2    44    0.0467  0.9544
plot(m1, main = "Residuals Shannon alpha diversity")

m2<-lme(Simpson~Drechslera_treatment*Fungicide_treatment*Cultivar,                            ,
        random=~1|Block,
        na.action=na.exclude,
        data=richness)

anova(m2)
##                                                   numDF denDF   F-value p-value
## (Intercept)                                           1    44 10405.612  <.0001
## Drechslera_treatment                                  1    44     1.177  0.2838
## Fungicide_treatment                                   1    44     0.129  0.7212
## Cultivar                                              2    44     0.631  0.5366
## Drechslera_treatment:Fungicide_treatment              1    44     0.070  0.7921
## Drechslera_treatment:Cultivar                         2    44     2.003  0.1471
## Fungicide_treatment:Cultivar                          2    44     1.206  0.3090
## Drechslera_treatment:Fungicide_treatment:Cultivar     2    44     1.526  0.2286
plot(m2, main = "Residuals Simpson alpha diversity")

Fungi

Loading data

setwd("C:/Users/yurip/Desktop/COMI/Colaborations/Stephan/ITS")

otu_mat<- as.data.frame(read_excel("asvs.ITS.xlsx"))
tax_mat<-  
  as.data.frame(read_excel("taxonomy.ITS.xlsx",
                           na = "NA",
                           col_names =
                             c("ASV", "Kingdom", "Phylum", "Class", "Order", "Family", "Genus", 
                                         "Species"))) %>%
  purrr::map_df(~gsub(
    pattern = "metagenome|uncultured|unidentified|Unknown",
    replacement = NA,
    .x)) %>%
  mutate_if(is.character, stringr::str_trim) %>%
  mutate(Kingdom = ifelse(is.na(Kingdom),
                          "U. Kingdom",
                          Kingdom),
         Phylum = coalesce(Phylum,
                           ifelse(grepl("^U.", Kingdom),
                                  Kingdom,
                                  paste("K_", Kingdom))),
         Class = coalesce(Class,
                          ifelse(grepl("^U.", Phylum),
                                 Phylum,
                                 paste("P_", Phylum))),
         Order = coalesce(Order,
                          ifelse(grepl("^U.", Class),
                                 Class,
                                 paste("C_", Class))),
         Family = coalesce(Family,
                           ifelse(grepl("^U.", Order),
                                  Order,
                                  paste("O_", Order))),
         Genus = coalesce(Genus,
                          ifelse(grepl("^U.", Family),
                                 Family,
                                 paste("F_", Family))),
         Species = coalesce(Species,
                            ifelse(grepl("^U.", Genus),
                                   Genus,
                                   paste("G_", Genus)))) %>%
  tibble::column_to_rownames("ASV") 

samples_df <- as.data.frame(read_excel("metadata_stephan.xlsx"))

otu_names=otu_mat$ASV
row.names(otu_mat) <- otu_names
otu_mat <- otu_mat[,2:61]


row.names(samples_df) <- samples_df$sample
samples_df <- samples_df %>% select (-sample) 
samples_df$Block <- factor(samples_df$Block, levels = c("Block1", "Block2", "Block3", "Block4","Block5" ))


otu_mat <- as.matrix(otu_mat)
tax_mat <- as.matrix(tax_mat)

class(otu_mat) <- "numeric"

OTU = otu_table(otu_mat, taxa_are_rows = TRUE)
TAX = tax_table(tax_mat)
samples = sample_data(samples_df)

stephan <- phyloseq(OTU, TAX, samples)
stephan
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1082 taxa and 60 samples ]
## sample_data() Sample Data:       [ 60 samples by 11 sample variables ]
## tax_table()   Taxonomy Table:    [ 1082 taxa by 7 taxonomic ranks ]
rank_names(stephan)
## [1] "Kingdom" "Phylum"  "Class"   "Order"   "Family"  "Genus"   "Species"
sample_variables(stephan)
##  [1] "Plant_number"             "Cultivar"                
##  [3] "Fungicide_treatment"      "Drechslera_treatment"    
##  [5] "Aphid_treatment"          "Block"                   
##  [7] "Cultivar_and_Drechslera"  "Fungicide_and_Drechslera"
##  [9] "Aphid_count"              "Total_Seed_Biomass"      
## [11] "Dry_Weight"

Rarefaction curve and normalization

mat <- t(otu_table(stephan))
class(mat) <- "matrix"
rarecurve(mat, step=50, cex=0.5, label = TRUE, 
          xlab = "Reads", ylab = "ASVs")

#Normalize number of reads in each sample using median sequencing depth.
stephan.2 <- subset_samples(stephan, !(Plant_number %in% c("27")))
total = median(sample_sums(stephan.2))
standf = function(x, t=total) round(t * (x / sum(x)))
stephan.rarefied = transform_sample_counts(stephan.2, standf)

Diversity analyzes

Alpha diversity

stephan.rarefied <- subset_samples(stephan.rarefied, !(Plant_number %in% c("27"))) ## I still need to put the new sample 27 in the dataset

alpha_fung=plot_richness(stephan.rarefied, x="Cultivar", color="Fungicide_treatment",
                        measures=c("Shannon")) + geom_boxplot() +
  facet_wrap(~Drechslera_treatment, scales = "free_x", nrow = 1) +
  scale_color_manual(values = c("#90EE90", "#228822")) +
  labs(title = "Alpha diversity (Shannon) - Fungi")

alpha_fung

alpha_fung2=plot_richness(stephan.rarefied, x="Cultivar", color="Fungicide_treatment",
                        measures=c("Simpson")) + geom_boxplot() +
  facet_wrap(~Drechslera_treatment, scales = "free_x", nrow = 1) +
  scale_color_manual(values = c("#90EE90", "#228822")) +
  labs(title = "Alpha diversity (Simpson) - Fungi")

alpha_fung2

## Calculation diversity metrics 

diversity=estimate_richness(stephan.rarefied)
shannon=diversity$Shannon
observed=diversity$Observed
observed.log <- log(observed)
evenness=shannon/observed.log
diversity$evennes=evenness
richness=as.data.frame(diversity)


DT::datatable(richness, class = "cell-border stripe", rownames = T,
              filter = "top", editable = TRUE, extensions = 'Buttons', options = list(
  dom = 'Bfrtip',
  buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
  scrollY = '300px',
  scrollCollapse = TRUE
))

Beta diversity

### PCoA and Constrained Analysis of Principal Coordinates (CAP)

ps1V4.rel_bray <- phyloseq::distance(stephan.rarefied, method = "bray") # CAP ordinate 
cap_ord <- ordinate(physeq = stephan.rarefied,  
                    method = "CAP", 
                    distance = ps1V4.rel_bray, 
                    formula = ~ Aphid_count + Total_Seed_Biomass +
                      Dry_Weight)

scree.cap <- plot_scree(cap_ord, "Scree Plot for MCs in Constrained Analysis of Principal Coordinates - Fungi")
print(scree.cap)

cap_plot <- plot_ordination(physeq = stephan.rarefied, 
                            ordination = cap_ord, 
                            color = "Cultivar", 
                            shape = "Fungicide_treatment",
                            axes = c(1,2)) + 
  geom_point(aes(colour = Cultivar), size = 3) + 
  geom_point(size = 3) +
  scale_color_manual(
    values = c( "#90EE90", "#228822", "#5F7FC7", "orange","#DA5724")) 



#Now add the environmental variables as arrows
cap_ord
## Call: capscale(formula = distance ~ Aphid_count + Total_Seed_Biomass +
## Dry_Weight, data = data)
## 
##               Inertia Proportion Rank
## Total         11.8754     1.0000     
## Constrained    1.3749     0.1158    3
## Unconstrained 12.0069     1.0111   34
## Imaginary     -1.5064    -0.1269   24
## Inertia is squared Bray distance 
## 
## Eigenvalues for constrained axes:
##   CAP1   CAP2   CAP3 
## 0.9305 0.3480 0.0964 
## 
## Eigenvalues for unconstrained axes:
##  MDS1  MDS2  MDS3  MDS4  MDS5  MDS6  MDS7  MDS8 
## 5.857 1.536 1.080 0.539 0.478 0.401 0.373 0.249 
## (Showing 8 of 34 unconstrained eigenvalues)
arrowmat <- vegan::scores(cap_ord, display = "bp")

# Add labels, make a data.frame
arrowdf <- data.frame(labels = rownames(arrowmat), arrowmat)

# Define the arrow aesthetic mapping
arrow_map <- aes(xend = CAP1, 
                 yend = CAP1, 
                 x = 0, 
                 y = 0, 
                 shape = NULL, 
                 color = NULL, 
                 label = labels)

label_map <- aes(x = 1.3 * CAP1, 
                 y = 1.3 * CAP1, 
                 shape = NULL, 
                 color = NULL, 
                 label = labels)

arrowhead = arrow(length = unit(0.02, "npc"))

# Make a new graphic
cap_plot + 
  geom_segment(
    mapping = arrow_map, 
    size = .5, 
    data = arrowdf, 
    color = "black", 
    arrow = arrowhead
  ) + 
  geom_text(
    mapping = label_map, 
    size = 4,  
    data = arrowdf, 
    show.legend = FALSE
  ) +
  stat_ellipse(type = "t") + ggtitle("CAP_Plot - Fungi")

Taxonomic composition

Box plot

mycols <- c("yellow", "#90EE90", "#228822", "#5F7FC7", "orange","#DA5724")

tax_table(stephan.rarefied)[,colnames(tax_table(stephan.rarefied))] <- gsub(tax_table(stephan.rarefied)[,colnames(tax_table(stephan.rarefied))],pattern="O_ C_ P_ K_ k__Fungi",replacement="Unidentified fungi") 


pn <- plot_taxa_boxplot(stephan.rarefied,
                        taxonomic.level = "Family",
                        top.otu = 9,
                        group = "Cultivar_and_Drechslera",
                        title = "top 9 taxa at family level - Fungi",
                        keep.other = FALSE,
                        add.violin= FALSE,
                        group.order = c("ACCO Drechs_NO", "ACCO Drechs_YES",
                                        "KLAR Drechs_NO", "KLAR Drechs_YES" ,
                                        "PLAN Drechs_YES", "PLAN Drechs_NO"),
                        group.colors = mycols
)
   
box_plot_fungi = pn + theme_biome_utils() +  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  theme(legend.position = "bottom")

box_plot_fungi

top_nine_fungi <- subset_taxa(stephan.rarefied, Family %in% c("Chaetomiaceae" ,  "Helotiaceae",   
                                                                 "Naviculisporaceae", "Nectriaceae" , 
                                                                 "Podosporaceae", "Unidentified fungi", "Pleurotheciaceae", "O_C_P_p__Ascomycota", "O_C_P_p__Chytridiomycota"))
  


top_nine_fungi
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 738 taxa and 59 samples ]
## sample_data() Sample Data:       [ 59 samples by 11 sample variables ]
## tax_table()   Taxonomy Table:    [ 738 taxa by 7 taxonomic ranks ]
top_nine_melted_fungi=psmelt(top_nine_fungi)


write_xlsx(top_nine_melted_fungi,"top_nine_fungi.xlsx")

Heatmap

meco_dataset <- phyloseq2meco(stephan.rarefied)

t1 <- trans_abund$new(dataset = meco_dataset, taxrank = "Family", ntaxa = 10, groupmean = "Fungicide_and_Drechslera")

t1 <- trans_abund$new(dataset = meco_dataset, taxrank = "Family", ntaxa = 20)
t1$plot_heatmap(facet = "Cultivar", xtext_keep = FALSE, withmargin = FALSE)

Statistical analyzes

Permanova

stephan_bray <- phyloseq::distance(stephan, method = "bray")
stephandf <- data.frame(sample_data(stephan))

adonis2(stephan_bray ~ Fungicide_treatment*Drechslera_treatment*Cultivar, data = stephandf)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = stephan_bray ~ Fungicide_treatment * Drechslera_treatment * Cultivar, data = stephandf)
##                                                   Df SumOfSqs      R2       F
## Fungicide_treatment                                1   0.1234 0.00982  1.0979
## Drechslera_treatment                               1   0.1537 0.01223  1.3679
## Cultivar                                           2   6.1188 0.48684 27.2221
## Fungicide_treatment:Drechslera_treatment           1   0.2027 0.01613  1.8035
## Fungicide_treatment:Cultivar                       2   0.2036 0.01620  0.9057
## Drechslera_treatment:Cultivar                      2   0.2158 0.01717  0.9600
## Fungicide_treatment:Drechslera_treatment:Cultivar  2   0.1559 0.01241  0.6937
## Residual                                          48   5.3946 0.42921        
## Total                                             59  12.5685 1.00000        
##                                                   Pr(>F)    
## Fungicide_treatment                                0.305    
## Drechslera_treatment                               0.199    
## Cultivar                                           0.001 ***
## Fungicide_treatment:Drechslera_treatment           0.126    
## Fungicide_treatment:Cultivar                       0.481    
## Drechslera_treatment:Cultivar                      0.436    
## Fungicide_treatment:Drechslera_treatment:Cultivar  0.682    
## Residual                                                    
## Total                                                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Linear mixed effect models

samples_df = samples_df %>% filter(Plant_number != 27)

richness$Block = samples_df$Block
richness$Drechslera_treatment = samples_df$Drechslera_treatment
richness$Fungicide_treatment = samples_df$Fungicide_treatment
richness$Sample = samples_df$Plant_number
richness$Cultivar = samples_df$Cultivar


table(richness$Fungicide_treatment,richness$Drechslera_treatment)
##                
##                 Drechs_NO Drechs_YES
##   Fungicide_NO         15         14
##   Fungicide_YES        15         15
m1<-lme(Shannon~Drechslera_treatment*Fungicide_treatment*Cultivar,                            ,
        random=~1|Block,
        na.action=na.exclude,
        data=richness)

anova(m1)
##                                                   numDF denDF   F-value p-value
## (Intercept)                                           1    43 1755.2338  <.0001
## Drechslera_treatment                                  1    43    0.0000  0.9961
## Fungicide_treatment                                   1    43    0.0005  0.9822
## Cultivar                                              2    43    0.2101  0.8113
## Drechslera_treatment:Fungicide_treatment              1    43    0.1546  0.6961
## Drechslera_treatment:Cultivar                         2    43    0.1601  0.8525
## Fungicide_treatment:Cultivar                          2    43    1.7716  0.1823
## Drechslera_treatment:Fungicide_treatment:Cultivar     2    43    1.1499  0.3262
plot(m1, main = "Residuals Shannon alpha diversity")

m2<-lme(Simpson~Drechslera_treatment*Fungicide_treatment*Cultivar,                            ,
        random=~1|Block,
        na.action=na.exclude,
        data=richness)

anova(m2)
##                                                   numDF denDF   F-value p-value
## (Intercept)                                           1    43 2252.2469  <.0001
## Drechslera_treatment                                  1    43    0.3411  0.5622
## Fungicide_treatment                                   1    43    0.0576  0.8115
## Cultivar                                              2    43    0.2144  0.8079
## Drechslera_treatment:Fungicide_treatment              1    43    0.0617  0.8050
## Drechslera_treatment:Cultivar                         2    43    0.2077  0.8132
## Fungicide_treatment:Cultivar                          2    43    1.2332  0.3015
## Drechslera_treatment:Fungicide_treatment:Cultivar     2    43    0.6009  0.5529
plot(m2, main = "Residuals Simpson alpha diversity")